using System;

using Atomic.Libraries.Mathematics.LinearAlgebra.AlgLib;

namespace Atomic.Libraries.Mathematics.LinearAlgebra
{
	[Serializable]
	public sealed class SingularValueDecomposition
	{
		private int m, n;

		private SingularValueDecomposition(Matrix a, bool full)
		{
			m = a.Rows;
			n = a.Columns;

			double[,] a0 = a.ToArray();
			double[] w = null;
			double[,] u = null;
			double[,] vt = null;

			int needed = full ? 2 : 1;
			svd.rmatrixsvd(a0, m, n, needed, needed, 2, ref w, ref u, ref vt);

			Full = full;
			A = a;
			U = new Matrix(u);
			V = Matrix.Transpose(new Matrix(vt));
			VTransposed = new Matrix(vt);
			W = new Vector(w);
			Sigma = Matrix.Diagonal(w);
		}
		
		public static SingularValueDecomposition Decompose(Matrix a)
		{
			return Decompose(a, false);
		}
				
		public static SingularValueDecomposition Decompose(Matrix a, bool full)
		{
			return new SingularValueDecomposition(a, full);
		}

		public static SingularValueDecomposition DecomposeFull(Matrix a)
		{
			return Decompose(a, true);
		}

		public Matrix PseudoInverse()
		{
			if (!Full)
			{
				throw new InvalidOperationException();
			}

			// http://en.wikipedia.org/wiki/Singular_value_decomposition#Pseudoinverse

			double[,] s = new double[n, m];
			for (int i = 0; i < m && i < n; i++)
			{
				double x = Sigma[i, i];
				if (x != 0.0)
				{
					s[i, i] = 1.0 / x;
				}
			}

			return V * new Matrix(s) * Matrix.Transpose(U);
		}

		public static Matrix PseudoInverse(Matrix a)
		{
			return DecomposeFull(a).PseudoInverse();
		}

		public double PseudoDeterminant()
		{
			if (m != n)
			{
				throw new ArgumentException("The pseudo-determinant is only defined for square matrices.");
			}

			// http://en.wikipedia.org/wiki/Pseudo-determinant
			// https://stat.ethz.ch/pipermail/r-help/2001-May/012909.html

			double d = 1.0;
			for (int i = 0; i < n; i++)
			{
				d *= Sigma[i, i];
			}

			return d;
		}

		public static double PseudoDeterminant(Matrix a)
		{
			return Decompose(a).PseudoDeterminant();
		}

		public bool Full
		{
			get;
			private set;
		}

		public Matrix A
		{
			get;
			private set;
		}

		public Matrix U
		{
			get;
			private set;
		}

		public Vector W
		{
			get;
			private set;
		}

		public Matrix Sigma
		{
			get;
			private set;
		}

		public Matrix VTransposed
		{
			get;
			private set;
		}
			
		public Matrix V
		{
			get;
			private set;
		}
	}
}
